En este tutorial vamos a ver cómo llevar a cabo un análisis de autocorrelación espacial usando tanto R como ArcGIS Pro y QGis.
Esto nos permitirá comparar las diferentes herramientas, analizar en qué difieren y en qué se parecen, y tendremos por tanto un mayor abanico de opciones para realizar este tipo de análisis.
En teoría hemos visto que la autocorrelación espacial no es más que el grado de correlación de una variable consigo misma a través del espacio. Para poder calcularla necesitamos por tanto: 1) observaciones de una variable, y 2) la ubicación espacial de dichas observaciones. Todo ello nos lo proporciona cualquier fichero de información espacial (shapefiles de ESRI, geopackages, u objetos de R de clase sf, como vimos en el tutorial “Accediendo a datos espaciales en R”).
Por ello, lo primero que vamos a hacer es cargar el fichero de datos espaciales que vamos a analizar, que no es otro que el shapefile con los valores de estaciones meteorológicas que usamos en la unidad sobre regresión lineal. También necesitamos cargar el paquete sfdep, que contienen funciones específicas para análisis espaciales en R. Por último, cargaremos el shapefile de provincias para mejorar nuestra visualización:
# Cargamos las librerías que necesitamoslibrary(sf)# library(spdep)library(sfdep) # Cargamos el fichero de datos (en este caso, un shapefile)estaciones <-st_read('data/meteo/meteo_espacial/estaciones_meteo.shp')
Reading layer `estaciones_meteo' from data source
`C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\estaciones_meteo.shp'
using driver `ESRI Shapefile'
Simple feature collection with 90 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N
Reading layer `provincias_spain' from data source
`C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\provincias_spain.shp'
using driver `ESRI Shapefile'
Simple feature collection with 51 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -14129.47 ymin: 3892590 xmax: 1126923 ymax: 4859517
Projected CRS: ETRS89 / UTM zone 30N
7.2 Definiendo vecinos
Hemos visto en la parte de teoría que uno de los pasos críticos en un análisis de autocorrelación es decidir qué observaciones consideramos como vecinas de un determinado punto. Hay varias manera de definirlo, pero las más importantes son definir los vecinos con criterios de contigüidad (en el caso de polígonos o raster), seleccionar los \(k\) puntos más cercanos, o definir el vecindario en función de la distancia entre los elementos. Métodos más complejos permiten también asignar pesos a los vecinos en función de diversos criterios, entre los que destaca asignarlos de manera inversa a la distancia entre puntos.
7.2.1 Vecinos contiguos
Si tenemos un fichero de polígonos, podemos definir la vecindad como aquellos polígonos contiguos a cada uno, mediante la función st:contiguity() del paquete sfdep. La opción queen = TRUE determina cómo vecinos de un polígono todos aquellos que coincidan en un sólo punto de sus límites. Si queen = FALSE sólo son vecinos aquellos que compartan al menos dos puntos (vecindad de Rook)
Neighbour list object:
Number of regions: 51
Number of nonzero links: 222
Percentage nonzero weights: 8.535179
Average number of links: 4.352941
4 regions with no links:
7 49 50 51
5 disjoint connected subgraphs
Ahora podemos plotear este objeto para ver las conexiones entre polígonos:
plot(st_geometry(provincias))plot(veins_cont, coords =st_coordinates(st_centroid(provincias)), add =TRUE, col ="red")
Warning: st_centroid assumes attributes are constant over geometries
NOTA: este método sólo funciona para archivos espaciales de polígonos. Si tratáramos de aplicarlos con puntos obtendríamos un error, ya que los puntos no pueden tener puntos de contacto entre ellos.
NOTA (2): aunque la opción “queen” permite definir el criterio de selección de vecinos, en el caso de polígonos complejos como el del fichero provincias lo más normal es que no haya diferencias entre queen = TRUE o queen = FALSE
7.2.2\(k\) vecinos más próximos
Un método alternativo consiste en elegir como vecinos los \(k\) puntos más cercanos, lo que se adapta a toda la zona de estudio, teniendo en cuenta las diferencias en las densidades de las entidades de área. En este caso tan sólo debemos definir el valor de \(k\), y la función buscará automáticamente los k puntos más próximos de cada punto. Usaremos para ello la función st_knn():
veins_k <-st_knn(estaciones, k =8)veins_k
Neighbour list object:
Number of regions: 90
Number of nonzero links: 720
Percentage nonzero weights: 8.888889
Average number of links: 8
Non-symmetric neighbours list
Ahora podemos plotearlo como antes:
plot(veins_k, coords =st_coordinates(estaciones), col ="red")
En el caso de esta metodología, puede utilizarse también con polígonos. Sin embargo requiere que trabajemos con los centroides de cada polígono, que serán los que nos permitan determinar la distancia entre polígonos:
veins_prov_k <-st_knn(st_centroid(provincias), k =4)
Warning: st_centroid assumes attributes are constant over geometries
Warning: st_centroid assumes attributes are constant over geometries
7.2.3 Vecinos basados en intervalos de distancia
Para definir los vecinos en base a una distancia podemos usar la función st_dist_band() en la que debemos definir la distancia mínima (lower) y máxima (upper) para seleccionar un vecino. La función determinará que todos los puntos que estén entre estas dos distancias serán vecinos:
veins_dist <-st_dist_band(estaciones, lower =0, upper =30000)plot(veins_dist, coords =st_coordinates(estaciones), col ="red")
Igual que antes, si queremos usar este método con polígonos, debemos hacerlo con los centroides de los polígonos.
7.3 Asignando pesos a los vecinos
Otra opción muy interesante, y que suele aplicarse en muchos casos, es no sólo definir cuáles consideraremos como vecinos de cada observación, sino darles un peso determinado según la distancia a la que estén del punto analizado. Este método se conoce como weighting.
Aunque existen varias opciones para dar pesos, veremos las dos más comunes: asignar el mismo peso a todos los vecinos de un punto (función st_weights()) o asignar mayor peso a los vecinos que estén más cerca, lo que se conoce como inverse distance weighting (idw) (función st_inverse_distance()).
Para evitar que la función calcule las distancias entre todos los puntos dos a dos - lo que podría colapsar la memoria del ordenador - las dos funciones se aplican sobre un objeto de vecinos ya creado, en el que le digamos cuáles son los vecinos de cada punto. En este caso, usaremos el que acabamos de crear (veins_dist):
Al imprimir el resultado vemos que, mientras en veins_dist_w() todos los vecinos de una observación tienen el mismo peso, en veins_idw los pesos varían entre los vecinos.
7.4 Autocorrelació Global: I de Moran
La I global de Moran determina el grado de autocorrelación de una muestra en su conjunto. Valores positivos indican autocorrelación positiva (los valores más próximos entre sí muestran valores más similares), valores negativos muestran lo contrario, y valores cercanos a 0, ausencia de autocorrelación espacial.
7.4.1 I global de Moran con R
En R podemos determinar la I de Moran con la función global_moran_test() que necesita como argumentos: (1) la variable que queremos analizar (en este caso será la temperatura media del mes, Tmed_MES); (2) un objeto neighbor que identifique los vecinos de cada observación; y (3) un objeto de pesos. Por lo tanto los pasos a seguir son:
Identificar los vecinos de cada observación
Asignar los pesos mediante la opción st_weights() o st_inverse_distance()
Calcular la I de Moran
Aunque ya tenemos los objetos resultantes de los pasos 1 y 2, vamos a repasar todo el proceso para tener todos los pasos juntos:
7.4.1.1 I global de Moran en base a los k vecinos más próximos
Como hemos visto antes, en el caso de usar los k vecinos más próximos, usaremos las funciones st_knn() y st_weights():
# 1. Identificar los vecinos de cada observación en un objeto neighbour (nb)veins_k <-st_knn(estaciones, k =8)# 2. Asignar los pesos de cada observación (en este caso, pesos iguales)veins_w_k <-st_weights(veins_k)# 3. Calcular la I de Moran para la variable Tmed_MESglobalMoran_k <-global_moran_test(x = estaciones$Tmed_MES, nb = veins_k, wt = veins_w_k)globalMoran_k
Moran I test under randomisation
data: x
weights: listw
Moran I statistic standard deviate = 12.934, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.607761749 -0.011235955 0.002290385
Vemos que el valor de la I de Moran es 0.6077, lo que indica autocorrelación espacial positiva. El valor esperado si los datos se distribuyeran aleatoriamente sería de -0.011. Por último, también sabemos que la probabilidad de que los valores de temperatura observados estén distribuidos al azar es de 2.2e-16, con lo que claramente rechazaremos la hipótesis nula (es decir, rechazamos que no haya autocorrelación). En definitiva, que sí que hay autocorrelación espacial.
7.4.1.2 I global de Moran en base a la distancia entre vecinos
En este caso reconoceremos como vecinos aquellos dentro de un intervalo de distancias, y asignaremos el mismo peso a todos:
# 1. Identificar los vecinos de cada observación en un objeto neighbour (nb)veins_dist <-st_dist_band(estaciones, lower =0, upper =30000)# 2. Asignar los pesos de cada observación (en este caso, pesos iguales)veins_w_dist <-st_weights(veins_dist)# 3. Calcular la I de MoranglobalMoran_dist <-global_moran_test(x = estaciones$Tmed_MES, nb = veins_dist, wt = veins_w_dist)globalMoran_dist
Moran I test under randomisation
data: x
weights: listw
Moran I statistic standard deviate = 11.134, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.664310960 -0.011235955 0.003681663
Vemos que el valor cambia respecto al calculado antes (0.664),ya que la definición de cuáles son los vecinos ha cambiado mucho. Sin embargo, la conclusión sigue siendo la misma (hay bastante autocorrelación)
7.4.1.3 I Global de Moran en base a IDW
Ahora identificaremos los vecinos en base a la franja de distancias, pero les asignaremos un peso inversamente proporcional a la distancia con el punto analizado:
# 1. Identificar los vecinos de cada observación en un objeto neighbour (nb)veins_dist <-st_dist_band(estaciones, lower =0, upper =30000)# 2. Asignar los pesos de cada observación (en este caso, pesos iguales)veins_idw <-st_inverse_distance(veins_dist, estaciones)# 3. Calcular la I de MoranglobalMoran_idw <-global_moran_test(x = estaciones$Tmed_MES, nb = veins_dist, wt = veins_idw)globalMoran_idw
Moran I test under randomisation
data: x
weights: listw
Moran I statistic standard deviate = 8.7647, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.663803406 -0.011235955 0.005931776
Aunque el resultado es parecido al de globalMoran_dist, hay pequeñas diferencias debidas a los pesos asignados. En función del dataset concreto las diferencias pueden ser más o menos grandes.
7.4.2 I Global de Mo#ran con ArcGis Pro
Determinar la I global de Moran con ArcGIS es también muy sencillo. En primer lugar, debemos cargar la capa estaciones_meteo.shp, y tenemos que buscar la función Spatial Autocorrelation (Global Moran’s I), que podemos encontrar dentro de las herramientas de Geoprocessing
Una vez seleccionada la función, tenemos que indicar la capa que queremos analizar, la variable de interés, y cómo determinaremos el peso de los vecinos.
Los resultados se muestran en la pestaña Results de ArcGIS, que si no tenemos visible, podemos activar yendo a Geoprocessing/Results. Si además hemos seleccionado que queremos que genere un informe, este estará disponible como html, produciendo una salida visualmente muy informativa.
> ** Ejercicio** Prueba a modificar la manera de definir los pesos de cada observación, verás que los resultados de ArcGIS son los mismos que obtenemos en R.
7.4.3 I Global de Moran con QGis y Geoda
7.4.3.1 Presentando GeoDa
Por supuesto, QGis tiene funcionalidades parecidas a las de ArcGis, mediante un plugin llamado Spatial Analysis Toolbox. Sin embargo, a día de hoy (julio de 2024), este plugin no funciona, así que veremos una alternativa muy interesante, el software GeoDa.
GeoDa es un software de escritorio gratuito y de código abierto, desarrollado por el equipo del Dr. Luc Anselin, profesor de la Universidad de Chicago y uno de los grandes pioneros en el análisis espacial y la econometría. Podemos descargar gratuitamente GeoDa aquí.
Una vez descargado, si lo abrimos nos pedirá que carguemos los datos: podemos cargar capas y objetos espaciales en multitud de formatos. En este caso, vamos a cargar el fichero estaciones_meteo.shp.
Una vez cargados los datos, el aspecto es el de un SIG convencional, con la tabla de contenidos a la izquierda, la visualización a la derecha, y arriba veremos el menú con diversas opciones.
7.4.3.2 Calculando vecinos y pesos en GeoDa
Para poder calcular la I de Moran, GeoDa nos exige que calculemos los pesos. Para ello debemos acceder al menú Weights Manager caracterizado por una gran W morada. El menú del Weights manager nos resultará familiar, ya que nos da la opción de definir los vecinos por contiguidad (según el criterio de Rooke o de Queen) o por distancia (y en este caso podremos elegir por banda de distancia, k vecinos más próximos o mediante un kernel, que no lo hemos visto). En caso de elegir los vecinos por distancia, nos preguntará si queremos usar el inverso de la distancia como peso (es decir, un idw):
Una vez decidido el criterio que queramos, guardamos el fichero de pesos en disco, con el nombre que queramos y la extensión .gwd. Una de las funcionalidades más interesantes de GeoDa es que nos permite visualizar, una vez seleccionado el fichero de pesos que queramos, los vecinos de un punto sobre el que pongamos el cursos, además de permitirnos ver gráficos de vecindad interactivos:
Podemos probar a crear un fichero con knn (8 vecinos), otro con los vecinos entre 0 y 30.000 metros, y otro con los mismos vecinos pero pesos inversos a la distancia, como hicimos en R.
7.4.3.3 Calculando la I de Moran con GeoDa
Una vez definidos los vecinos, podemos calcular fácilmente la I global de Moran, seleccionando el comando Univariate Moran's I del menú Space. Para ello sólo debemos indicar la variable de interés (Tmed_MES) y el fichero de vecinos deseado:
7.5 Autocorrelación Local (I de Anselin)
Para poder calcular la I local de Moran (también llamada I de Anselin), necesitaremos, igual que antes, conocer la variable a utilizar, la identificación de los vecinos, y los pesos. Vamos a trabajar, por simplicidad, con uno de los casos (idw), pero funcionaría igual con el resto:
localMoran <-local_moran(x = estaciones$Tmed_MES, nb = veins_dist, wt = veins_idw)
Vemos que en este caso el objeto generado contiene el valor calculado de I (ii), el valor esperado (eii), la varianza (var_ii), el z-score (z_ii) y el p-valor (p_ii) de cada observación de la muestra, entre otros campos
7.5.1 LISA (Local indicators of spatial association) - Agrupaciones espaciales
Un análisis muy interesante cuando trabajamos con autocorrelación espacial local es el de identificar aquellas observaciones en las que se cumplen una serie de condiciones. Este análisis fue desarrollado en 1995 por Luc Anselin, y recibe el nombre de LISA (local indicators of spatial association). Se compone de varios elementos:
7.5.1.1 Scatterplot de Anselin
Si queremos producir un scatterplot con los valores positivamente y negativamente autocorrelacionados, podemos usar la función moran.plot(), del paquete spdep. Por desgracia, aún no existe una versión de esta función para sfdep, por lo que nos vemos obligados a definir los pesos de los vecinos usando la función nb2listw() (en caso de querer aplicar pesos iguales), o la función nb2listwdist() si queremos aplicar un idw.
library(spdep)
Cargando paquete requerido: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Vemos que los valores de las esquinas inferior izquierda y superior derecha son los que tienen autocorrelación positiva, y los de las esquinas contrarias, autocorrelación negativa.
7.5.1.2 Mapa de valores de I
También podemos usar la función tm_shape() y tm_dots() del paquete tmap para visualizar espacialmente los valores locales de I de Moran. Primero debemos unir los cálculos de I local al objeto espacial que tenemos, usando cbind().
lisa <-cbind(estaciones, localMoran)
Después podemos crear la representación espacial de la siguiente manera:
library(tmap)
Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
Variable(s) "ii" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
7.5.1.3 Mapa de agrupaciones cluster
Consiste en identificar las observaciones que cumplen estos criterios:
Alto valor de la variable y vecinos con alto valor (High-high)
Bajo valor de la variable y vecinos con bajo valor (Low-Low)
Alto valor de la variable y vecinos con bajo valor (High-Low)
Bajo valor de la variable y vecinos con alto valor (Low-High)
Los dos primeros se suelen colorear de color rojo y azul intenso, respectivamente, y corresponden a las observaciones que contribuyen significativamente a la autocorrelación espacial global positiva. Los otros dos se suelen representar con colores pálidos, y son las observaciones que contribuyen a una autocorrelación negativa. Las observaciones con valores no significativos no suelen colorearse.
La función localMoran ya nos ofrece una clasificación en estos cuatro grupos, a la que podemos acceder de la siguiente manera:
Como antes hemos unido los resultados de localMoran a nuestra tabla estaciones también tendremos esta información en el objeto moran.map
Si definimos un nivel de significación, podemos clasificar las observaciones según el grupo al que pertenezcan y su significación. Indiquemos que los valores no significativos tengan un valor de “NS”
# llindar de significacionsignif <-0.1lisa$mean <-ifelse(lisa$p_ii > signif, "NS", as.character(lisa$mean))lisa
Por último, ploteamos la figura con tmap, tm_shape() y tm_dots():
# plot dels resultatstm_shape(lisa) +tm_dots(col ="mean",style ="cat",palette=c("NS"="grey","Low-Low"="blue", # Low-Low"Low-High"="lightblue", # Low-High "High-Low"="sienna1", # High-low"High-High"="red"), # High -Highsize =0.5) +tm_shape(provincias)+tm_borders()
7.5.1.4 Mapa de significación
Otro elemento interesante es producir un mapa de significación. Para ello establecemos las categorías de significación que queramos, y ploteamos su distribución espacial. Por ejemplo, vamos a crear 4 categorías:
> ** Ejercicio** Prueba a modificar la manera de definir los pesos de cada observación, verás que los resultados de ArcGIS son los mismos que obtenemos en R.
Una vez decidido el criterio que queramos, guardamos el fichero de pesos en disco, con el nombre que queramos y la extensión 
